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NOMENCLATURE 


Text 

A  coefficient  matrix 

al...a5  coefficients 

cO  coefficient  matrix  of  dependent  variable  in  the  finite-volume  discretization 

cl  coefficient  matrix  in  x  direction  of  the  differential  term  in  finite-volume  equation 

c2  coefficient  matrix  in  y  direction  of  the  differential  term  in  finite-volume  equation 

E  total  energy  density  of  gas  mixture,  ergs/cm^ 

f  source  term  vector 

H  height,  cm 

J  molecular  diffusive  flux,  gmole/cm^sec 

L  length,  cm 

P  pressure,  dynes/cm^ 

6P  change  in  pressure,  dynes/cm^ 

q  conductive  heat  flux,  ergs/cm^sec 

r,s  residual  vectors 

T  temperature,  K 

t  time,  sec. 

u  dependent  variable  scaler 

u  dependent  variable  vector 

V  velocity  vector 

W  rate  of  combustion,  gmole/cm^sec 

X  horizontal  axial  coordinate 

y  vertical  axial  coordinate 

Greek  Symbols 
p  gas  density,  kg/m^ 

8  energy  density  rate,  ergs/cm^-sec. 

Y  specific  heat  ratio 

T  stress  tensor 

AHc  heat  of  combustion 

to  implicitness  parameter 

Subscripts 

i  index  in  x-direction 

j  index  in  y-direction 

k  k  th  specie 

r  radiative  component 

Superscripts 
n  new  time  step 

o  previous  time  step 


NUMERICAL  SIMULATION  OF  SOLID  COMBUSTION 
WITH  A  ROBUST  CONJUGATE-GRADIENT  SOLUTION  FOR  PRESSURE 


1.  INTRODUCTION 

Solid  materials,  such  as  the  insulation  on  bulkheads  and  doors  in  a  ship,  often  form 
boundary  layer  flames  due  to  no-slip  and  pyrolysis  of  the  solid.  The  burning  rate  plays  a  key 
role  in  the  evaluation  of  fire  threat  and  depends  on  both  the  fluid  dynamics  near  the  surface  and 
the  intrinsic  burning  characteristics  of  the  solid  phase.  A  numerical  model  was  developed  by 
Ananth  et  al,  [1-4]  to  describe  temperature  and  local  burning  rate  distributions  along  the  length 
of  a  polymer  (poly  methyl  methacrylate,  PMMA)  slab.  Time  dependent  solutions  of  full 
Navier-Stokes  equations  are  obtained  by  using  Barely  Implicit  Correction  to  Flux  Corrected 
Transport  (BIC-FCT)  algorithms  for  combustion  of  a  PMMA  plate  under  forced  flow  of  air  past 
the  plate.  A  multi-variable  fixed  point  (MVFP)  iterative  method  is  developed  to  describe  the 
coupling  between  the  gas  and  the  solid  phases. 

to  BIC-FCT  schemes,  an  implicit  elliptic  partial  differential  equation  for  pressure  arises 
from  the  momentum  and  energy  equations  as  described  by  G.  Patnaik  et  al.  [5].  In  the  simulation 
of  solid  combustion,  the  elliptic  equation  must  be  solved  at  every  MVFP  iteration  and  at  each 
time  step.  Since  the  solution  of  the  pressure  equation  is  of  significant  computational  cost,  a 
robust  method  is  needed.  The  objective  of  this  work  is  to  implement  a  bi-conjugate  gradient 
method  for  the  solution  of  the  pressure  equation  and  predict  the  burning  rate  distributions  along 
the  polymer  slab.  We  will  show  that  a  bi-conjugate  gradient  method  (Bi-CGSTAB)  is 
computationally  reliable,  efficient,  and  accurate  for  slow  flow  past  the  solid  plate  by  comparing 
it  with  a  multigrid  method  (MGRID),  which  is  used  currently.  Unlike  MGRID,  Bi-CGSTAB  is 
simple,  transparent,  easier  to  implement  with  direct  accuracy  control,  and  can  port  among 
different  machines  (e.g.,  SGI  unix  work  station  and  the  SunSPARC  cluster). 

To  understand  how  the  Bi-CGSTAB  method  works,  it  is  tested  using  a  simple  problem, 
which  has  an  analytical  solution  prior  to  the  simulation  of  the  combustion  problem.  Both 
MGRID  and  Bi-CGSTAB  methods  are  applied  to  the  solid  combustion  problem,  and  the 
solutions  and  CPU  times  compared  at  steady  state  (12,000  time  steps).  The  method  of  Bi- 
CGSTAB  shows  favorable  features  in  convergence  and  accuracy. 


2.  LITERATURE  REVIEW 

Numerical  modeling  of  combustion  problems  often  involves  the  solution  of  elliptic 
partial  differential  equations  subjected  to  complicated  boundaries  and  interfaces.  The  general 
form  of  an  elliptic  equation  is  given  by 

ai  u  +  a2  V*  a3Vu  + Va4U  =  a5  (2.1) 
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with  boundary  conditions  given  by 

bi  u  +  ba  ^  =  ba .  (2.2) 

dn 

Here  u  is  a  dependent  variable,  ai  to  as  and  bi  to  bs  are  coefficients,  and  n  representing  normal  to 
the  boundary  can  be  x  or  y. 

In  general,  the  elliptic  partial  differential  equation  is  discretized  by  a  finite  volume 
method.  This  results  in  a  set  of  linear  algebraic  equations,  Au=f,  where  A  is  the  coefficient 
matrix  and  u  and  f  are  the  vectors  containing  the  dependent  variables  and  source  terms, 
respectively.  Since  the  elliptic  partial  differential  equation  is  linear,  it  leads  to  a  linear  set  of 
algebraic  equations,  which  may  be  solved  using  an  iterative  algorithm.  However,  in  time 
dependent  formulations,  such  as  BIC-FCT,  the  coefficients  in  the  elliptic  equation  are  time- 
dependent.  Therefore,  the  coefficients  contained  in  the  matrix  A  vary  with  time.  As  combustion 
progresses  in  time,  the  coefficient  matrix  may  rapidly  fill  up  and  become  non-symmetric  or  non- 
triangular.  In  these  situations,  common  elliptic  solvers,  such  as  a  direct  method,  and  the  iterative 
technique  ADI  (Alternating-Direction  Implicit)  method,  become  expensive  [6].  In  recent  years, 
multigrid  methods  and  conjugate  gradient  methods  have  been  recognized  as  feasible  methods  to 
solve  such  elliptic  problems  to  achieve  the  required  accuracy  and  efficiency. 

The  simplest  of  the  conjugate  gradient  methods  is  based  on  the  idea  of  minimizing  the 
function 


S(u)  =(1/2)  u»(Au)  -  f  •  u  (2.3) 

The  function  S  is  minimized  when  the  gradient  is  zero  as  shown  below 
VS  =  Au-f=0  (2.4) 

The  iterative  method  uses  successive  approximations  to  reach  the  final  solution.  Each 
successive  iteration  searches  for  a  correction  to  the  previous  iteration  in  a  direction  normal  to  the 
subspace  of  the  most  recent  approximation,  until  the  entire  vector  space  is  searched.  The 
methods  work  well  for  those  slow  flows  problems,  wherein  viscous  or  diffusive  effects  are 
dominant.  A  typical  conjugate  gradients  method  for  solving  elliptic  equations  is  Conjugate 
Gradients-Squared  method  (CG-S)  of  Sonneveld,  et  al.[7].  In  this  method,  the  residual  vectors  r 
=Au-f  generated  by  the  method  satisfy  a  3-term  recurrence  relation.  The  method  has  the  feature 
of  fast  convergence. 

Bi-CG  methods  are  developed  based  on  CG-S  and  do  not  directly  relate  to  the  fiinction 
minimization  concept.  These  methods  have  two  3 -term  recurrence  relations  for  better  stability. 

It  is  found  that  CG-S  based  methods  may  lead  to  a  rather  irregular  convergence  behavior  and  in 
some  cases  rounding  errors  can  even  result  in  severe  cancellation  effects  in  the  solution.  Van 
Der  Vorst  [8]  proposed  a  variant  of  Bi-CG,  named  Bi-CGSTAB,  to  eliminate  these  negative 
effects.  Van  Der  Vorst  showed  that  the  Bi-CGSTAB  provides  fast  and  smoothly  convergence 
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without  giving  up  the  attractive  speed  of  convergence  of  CG-S.  G.  Patnaik  [9]  has  recently 
implemented  the  Bi-CGSTAB  method  to  study  3-dimensional  premixed  hydrogen-air  flames. 

He  also  successfully  parallelized  the  code  to  run  on  0rigin2000  machine.  The  present  work 
implements  the  Bi-CGSTAB  to  the  solid  combustion  problem  and  is  an  extension  of  the  work  of 
G.  Patnaik  [9]. 

We  have  been  using  a  multigrid  method  (MGRID)  developed  by  DeVore  [10]  to  solve 
the  pressure  equation.  This  method  for  solving  the  elliptic  boundary-value  problem  was 
originally  found  by  Douglas  [11].  Subsequently,  DeVore  [10]  converted  the  method  to 
vectorized  form  for  higher  speed  computations.  Basically,  the  method  uses  the  solution  obtained 
from  a  sequence  of  smaller  spaces  to  approximate  the  desired  solution  in  the  largest  space.  Since 
the  smaller  spaces  have  geometrically  fewer  unknowns,  they  require  less  computational  effort  to 
yield  a  solution  than  does  the  largest  space.  The  method  combines  direct  and  iterative  methods. 
(It  can  make  use  of  coarse  grids  with  mesh  spacing  for  two  fine-grid  points,  say  2h  as  two  levels, 
4h,  8h,  or  more  levels.)  The  approximate  solution  for  the  fine-grid  can  be  obtained  by  a  direct 
method,  such  as  Gaussian  elimination.  In  order  to  improve  the  approximated  solution  on  the  fine 
grids,  the  method  does  a  small  number  of  relaxation  iterations  on  the  solution  obtained  from  the 
first  approximation.  This  iteration  process  is  based  on  the  residual  vector,  and  then  the  desired 
values  from  the  fine  grid  are  projected  onto  the  coarse  grid.  This  process  is  a  correction  process, 
or  a  cycle,  and  can  be  done  by  several  iteration  steps.  According  to  the  report  from  DeVore  [10], 
a  level-4  and  cycle-2  algorithm  might  be  a  good  choice.  However,  he  recommended  that  the  best 
choices  for  many  of  these  options,  in  general,  depend  on  the  properties  of  the  coefficients  in  the 
equation,  the  number  of  spatial  dimensions,  the  number  of  unknowns,  and  the  accuracy  required 
in  the  solution.  Thus,  the  method  is  complex  and  less  transparent.  Furthermore,  for  the  solid 
combustion  problem,  the  results  reported  here  will  show  that  it  is  also  time  consuming,  due  to  the 
intergrid  transfers  and  recursive  cycling. 


3.  MATHEMATICAL  FORMULATION 

We  consider  air  flow  past  a  PMMA  plate,  which  includes  a  non-burning  leading  plate  and 
a  post  flame  plate  as  shown  in  Figure  1.  The  plate  is  ignited  near  the  entire  surface  by  adding 
external  energy  for  a  short  time.  The  fuel  vapor  generated  at  the  surface  mixes  with  the 
incoming  air  by  convection  and  diffusion  and  undergoes  combustion.  The  development  of  a 
boundary  layer  diffusion  flame  is  described  by  the  time-dependent  mass,  momentum,  total 
energy,  and  specie  equations,  which  are  given  by 

continuity 

^  =  -V«pv,  (3.1) 

dt 
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Figure  1  Boundary  Layer  Combustion  of  a  Solid  Fuel  Plate 


momentum 


d  ov  =  -V«pvv  -  VP+V  •  X ,  (3.2) 

at 


energy 


as  =  -V*(E  +  P)  v+V  •  (v«  T )+  V  •  q+WkAHc+qr ,  (3.3) 

at 

and  specie 

^  -  V  •  vCk+V  •  Jk+Wk  ,  (3.4) 

at 

where  the  total  energy  density  E  is 
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(3.5) 


E  =  s  + l_p 
2 

The  pressure  P  and  internal  energy  density  s  are  related  by  the  equation  of  state,  which  is  given 
by 


P  =  (Y  -  1)  s-  (3.6) 

Here  p  is  the  gas  density,  v  is  the  velocity  vector,t  is  the  stress  tensor,  q  is  conductive  heat  flux, 
Wk  is  the  rate  of  combustion  of  specie  k,  AHc  is  the  heat  of  combustion,  qr  is  the  radiative  heat 
flux,  Ckis  the  specie  concentration,  Jk  is  the  molecular  diffusive  flux,  and  y  is  the  ratio  of  specific 
heats. 


The  convective,  diffusive,  and  reaction  parts  of  the  equations  (3.1)  to  (3.6)  are  computed 
separately  and  are  combined  using  time  step-splitting  [12].  The  pressure  term  is  included  in  the 
convection  algorithms,  which  are  subject  to  the  Courant  condition  in  an  explicit  formulation. 
The  convective  part  of  the  equations  (3.2)  and  (3.3)  are  given  by 

5pv  =  -V*pvv-VP.  (3.7) 

3t 


and 


^  =  -V«(E  +  P)v,  (3.8) 

at 

To  solve  this  time-dependent  system  of  convection  equations  numerically,  Casulli  and 
Greenspan  [13]  showed  that  it  is  not  necessary  to  treat  every  term  in  a  finite-difference  algorithm 
implicitly  to  avoid  the  time  step  constraint  imposed  by  the  Courant  condition.  They  showed  that 
only  those  terms  containing  the  pressure  in  equation  (3.7)  and  the  velocity  in  equation  (3.8)  must 
be  treated  implicitly.  Patnaik  et  al.  [5]  applied  the  concepts  of  Casulli  and  Greenspan  [13]  to 
obtain  solutions  for  pressure  and  the  velocity  in  by  BIC-FCT  formulation.  An  implicitness 
parameter,  co,  is  introduced  and  it  can  be  vary  the  degree  of  implicitness  of  the  algorithm.  In 
general,  we  can  have  0.5  <  ©  <  1,  where  the  implicit  terms  are  centered  in  time  for  ©  =  0.5.  For 
©  <  0.5,  the  method  is  foimd  to  be  unstable  for  sufficiently  large  time  steps. 

The  algorithm  has  two  stages.  One  stage  is  an  explicit  predictor  that  determines  the  provisional 
value  p  and  v  in  the  following  two  equations, 


D  -  0°  =  -  V  •  P°  V° 

At 


(3.9) 


pv  -p°v°  =  -v,p°vOvO_vp°.  (3  10) 

At 
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The  bar  denotes  predictor  values  at  the  new  time  step,  and  the  superscript  o  and  n  are  used  to 
denote  the  values  at  old  time  step  and  new  time,  respectively,  in  the  correction.  The  implicit 
forms  of  momentum  and  energy  equations  are 

p"v"  -p°vO:=-V,pOvOvO_V,  [eoP"  +(1-Q))P°]  (3.11) 

At 


E"  -E°  =  -V>rE°  +  P°H^v"  +(l-co)v°].  (3.12) 

At 

Since  co  =1,  the  algorithm  is  completely  implicit  and  reverts  to  the  original  equations  as 
equations  (3.7)  and  (3.8). 

The  change  in  pressure,  5P,  is  defined  as 

8P  =a)(P"-P°).  (3.13) 

Then  the  correction  equation  for  momentum  can  be  obtained  in  terms  of  8P  by  subtracting 
equation  (3.10)  from  (3.1 1), 

d"  v"  -  o  V  =  -  Vm  r  P”  -  P°  '1  =  -  V  8P  ■  (3.14) 

At 

By  rearranging  equation  (  3.14)  and  letting  p"  =  p,  the  new  velocity  is  obtained  as 

v"  =  -4t  V5P+  V  ,  (3.15) 

P 

then  a  correction  equation  for  energy  using  the  equation  of  state  is  obtained  as 

s"=  5P  +s°  ,  (3.16) 

(Y  - 1)  CO 

where  the  co  factor  appears  from  the  definition  of  6P  from  equation  (3.13).  5P  can  be  found  by 
substituting  equations  (3.15),  (3.16),  and  (3.5)  into  equation  (3.12)  to  obtain 

p  V 2  -  qO  v°2  +  SP  =  Oi  At  V#  (  E°  +  P°  )  V  6P 
2  At  (Y“1)®  p 

-  CO  V*(  +  P")  V  -  (1-©)  V*(  E°  +  P°)  v”  .  (3.17) 
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Note  that  Ae  kinetic  energy  change  is  included  explicitly.  For  convenience,  we  define  the 
quantity  E, 

_l^  =  _V.(E°  +  P°)[co  v  +  (l-®)  v°]  (3.18) 

At 

so  that  equation  (3.17)  can  be  rewritten  as 

8JP _ -  (0  At  V*  (  E  jf  P°  1 V 5P  =  1-E”-  o  -o°v°^ 

-(y-l)(oAt  p  At  2At  .  (3.19) 

By  dividing  -  ©  At  from  both  sides  of  Eq.  (3.19),  then 

5P  +V«(  E°  +  PMV5P=  E-E°+  pv^ 

-(y-l)®^At^  p  -®At^  2©At^  (3.20) 

The  barely  implicit  correction  (BIC)  is  carried  out  in  three  stages.  In  the  first  stage, 
equations  (3.9),  (3.10),  and  (3.18)  are  integrated  with  any  one-step  explicit  method.  In  the 
second  stage,  the  pressure  correction  equation  (3.20)  is  solved  by  a  numerical  elliptic  solver.  In 
the  last  stage,  equations  (3.15)  and  (3.16)  are  used  to  obtain  the  values  of  momentum  and  energy 
at  the  new  time-step.  These  values  of  the  momentum  and  energy  are  then  added  to  the 
contributions  from  the  diffusive,  reactive,  and  radiative  parts  of  the  equations  (3.2)  and  (3.3)  to 
obtain  final  values  of  the  energy  and  momentum. 

Our  focus  here  is  to  obtain  the  solution  of  equation  (3.20)  for  8P,  which  is  of  the  form 
given  by  equation  (2. 1).  The  right  hand  side  of  equation  (3.20)  is  evaluated  explicitly. 

Therefore,  the  right  hand  side  of  equation  (3.20)  is  known  and  the  only  unknown  is  5P.  Clearly 
equation  (3.20)  is  linear  and  can  be  solved  by  Bi-CGSTAB  iterative  method. 


4.  NUMERICAL  FORMULATION 

In  order  to  solve  5P  numerically,  equation  (3.20)  can  be  written  in  the  form  of  equation 
(2.3)  by  recognizing 


u=5P 


(4.1) 


where 


ai  =  -1 _ 

(y-1)©^  At^ 


(4.2) 


32  =  1  , 


(4.3) 
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(4.4) 


a3=  E°  +  P° 

P 

a4=0 ,  (4.5) 

and  as  -  E-E°  +  p  v  ^  -  o°  (4.6) 

-  (0  At  ^  2  toAt^ 

Therefore,  ai..a5  are  known  and  5P  is  the  unknown.  In  general,  at  a  fixed  time,  the  coefficients 
ai  ...a5  and  variable  5P  vary  with  position  (x  and  y).  Next,  a  set  of  algebraic  equations  are 
derived  from  equations  (4.1)  to  (4.6)  following  the  procedure  outlined  by  G.  Patnaik  [9]  for 
premixed  flames. 

4.1  Algebraic  Equations  in  Matrix  Form 

Equation  (4.1)  is  discretized  based  on  a  two  dimensional  finite  volume  grid.  A  typical 
cell  is  shown  in  Figure  2.  Here,  the  grid  contains  nx  and  ny  cells  in  x  and  y  directions, 
respectively.  Also,  i  and  j  are  cell  numbers  in  x  and  y  directions,  respectively. 

The  descritized  equation  for  a  typical  cell  (i  J)  is  given  by 

Co(i,j)u(i,j)+  [ci(ij)u(i-l  j)  -Ci(i+l,j)u(ij)]  +[ci(i+l,j)u(i+l,j)-  Ci(i,j)u(i,j)] 

+  [C2(i,j)u(ij-1)  -C2(i,j+l)u(i,j)]  +  [c2(i,j+l)u(i,j+l)-C2(i,j)u(i,j)  ]=  f(i,j) 

for  i=l,nx  and  j=l,ny  (4.7) 

The  values  of  the  coefficients  and  the  variables  at  cell  numbers  0,  (nx+1),  and  (ny+1) 
correspond  to  the  guard  cells  and  are  evaluated  fi-om  the  boundary  conditions.  Inside  the  domain 
they  are  given  by  the  following  equations.  The  first  term  (self-term)  coefficients  of  equation 
(4.7)  are  given  by 

Co(i,j)=Ax(i)AyG)ai(i,j).  (4.8) 

The  forcing  terms  on  the  right  hand  side  of  equation  (4.7)  are  given  by 
f(i,j)  =Ax(i)AyG)a5(i  j).  (4.9) 
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Figure  2  A  Typical  Cell  for  Bi-CGSTAB 

The  X  coefficients  in  equation  (4.7)  are  given  by 

ci(i  j)=  0.5(a3(i- 1  j)+a3(i  j))2AyO‘)/(  Ax(i- 1 )  )+Ax(i)).  (4.10) 

The  y  coefficientsin  equation  (4.7)  are  given  by 

C2(iJ)=  0.5(a3(iJ-l)+a3(i,j))2Ax(i)/(  Ay(j-1)  )+Ay(j)).  (4. 1 1) 

If  the  coefficient  matrix  a2  is  not  an  identity  matrix,  then  Ci(i  j)  and  C2(i  j)  become 

ci(i j)  =  a2(i j)ci(i j),  .  (4.12) 

and 

C2(iJ)  =  a2(i,j)c2(i,j).  (4.13) 

Therefore,  equation  (4.7)  can  be  rewritten  as 

[K(i,j)]u(iJ)  +  C](i+l,j)u(i+l  j)+  ci(i,j)u(i-l,j)  +C2(i,j)u(i,j-1)  +C2(i,j+l)u(i,j+l)  =  f(ij) 
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for  i=  1  ,nx  and  j= 1  ,ny, 


(4.14) 


where 

K=[co(ij)-Ci(iJ)-c,(i+lJ)-C2(i,j)-C2(iJ+l)] .  (4.15) 

Equations  (4.14)  and  (4.15)  represent  a  set  of  linear  algebraic  equations  of  the  form 

(4.16) 


which  are  solved  by  an  iterative  method  as  discussed  in  section  4.2. 

4.2  Iterative  Solution 


The  Bi-CGSTAB  algorithms  developed  by  Van  Der  Vorst  [8]  are  shown  in  Table  1 . 
Initially,  an  arbitrary  vector  r  is  determined,  say  by  r  =  ro.  where  ro  is  an  initial  residual  matrix. 
The  values  of  r©  are  calculated  from  r=f-Auo.  Then  the  iteration  loop  process  starts  (I=l,2,3...). 
Within  each  loop,  Bi-CGSTAB  carries  out  two  updates  to  the  residual  r  of  the  current  solution  u. 
The  two  updated  residuals  are  calculated  from  s  =  r,.]  -  avi  and  n  =s-o)i  t.  This  explains  Bi-  as 
the  name  of  the  algorithm.  If  the  residual  s  is  small  enough,  the  second  updated  residual  will  be 
skipped.  The  iteration  process  is  continued  until  a  pre-specified  tolerance  condition  is  satisfied. 

Preconditioning  or  reformulation  of  the  equation  (4.14)  is  often  useful  in  speeding  up  the 
convergence.  Van  Der  Vorst  suggested  the  following  preconditioning  of  the  equations 


r  = 

{u(i,j)  +  K-‘[  cl(i+l,j)u(i+l,j)+  cl(i,j)u(i-l,j)  +c2(i,j)u(i,j-l)  +c2(i,j+l)u(i,j+l)]}-K'’  f(ij). 


where 

K-'=l/[c0(i,j)-cl(i,j)-cl(i+l,j)-c2(i,j)-c2(ij+l)]. 


(4.17) 

(4.18) 


If  the  magnitude  of  K  is  less  than  the  magnitude  of  K’‘,  then  the  magnitude  of  residual 
vector  is  reduced  and  the  iteration  process  is  faster,  because  the  initial  guess  of  u  is  closer  to  the 
solution.  The  iteration  process  with  a  preconditioned  algorithm  will  not  change  the  original 
system.  However,  it  changes  r,  where  r  is  an  arbitrary  vector,  e  g.,  r  =to,  and  ro^f-Axo  for  initial 
guess. 


The  convergence  criterion  in  Bi-CGSTAB  is  determined  by  the  square  root  of  the 
residual,  i.e.,  a  =  square  root  of  (ro^  rj.]),  where  (ro^t  rj.i)  is  a  dot  product  of  two  residual  vectors 
from  sequential  steps  in  the  calculation.  The  tolerance  is  selected  based  on  this  convergence 
criterion  and  is  specified  directly.  In  the  computer  code,  a  parameter  named  norm  is  introduced 
outside  iteration  loop,  and  norm  =  1 /square  root  of  (f,f). 
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BKGSIABAIgoiifliii 

Initial  guess  u  ^>ib=f-AU) 

Initial  guess  U)  =^>ib =f-Aq) 

Ai  aibitraty  vectorr  =io 

Ai  artihaiy  vector  r  =ib 

pc»=a=Q)=l 

po=a=ai)=l 

Vo=P)=0 

Vo=P)=0 

Rr  i=  .  (start  iteratoi  tocps) 

RrF=l,%... 

P=feri-i) 

|3=(pi/pi.i)(a'ca.i) 

P=(p/pi-i)(a'ca.]) 

prri-i+B  fe-i-oj-iV-i) 

Fpri-i+B  ^1-ui^iVi.i) 

SdveyfioTiI^pp 

vfAp 

vfAy  (>^’) 

0FPi/fe),y) 

s=Ti.i  -cR^  (ipdateresidiail) 

S=^.l-0(Vi 

SdveyfcmK^ 

bAs 

bAz  (z^bK*) 

a)=(fsy(1,t) 

G)=ac't,K'sy(K't,K:'t) 

dr^U-i+cpr+QS 

U=Ui.i+Cp-4tQS 

ffu  is  accurate  enough  then  quit 

ffu  is  accurate  enou^  then  quit 

i;  =s-cQ  t  (ipdateiesiduah 

ii=s-c}t 

end 

end 

Table  1  Bi-CGSTAB  Algorithms  with  and  without  Pre-conditioning  [8] 


The  convergence  criterion  is  set  as  the  product  of  a  and  norm.  This  could  help  set  a 
suitable  criterion  to  avoid  unrealistic  tolerance,  since  usually  a  is  a  very  small  number  if  the 
program  converges.  However,  if  the  source  term  f  is  zero,  norm  cannot  be  used. 

As  discussed  earlier,  the  values  at  the  guard  cells,  which  are  on  the  boundaries  of  the 
domain,  appear  explicitly  in  the  governing  algebraic  equations  (4.14).  In  Bi-CGSTAB  method, 
the  boundary  conditions  are  specified  at  each  iteration  step  when  the  updated  residual  is 
calculated.  For  the  solid  combustion  problem,  the  pressure  gradient  is  set  to  zero  at  all 
boundaries  except  the  outflow.  At  the  downstream  of  the  fuel  plate,  we  set  pressure  to  be 
atmospheric,  i.e.,  5P=0.  In  BIC-FCT  formulation,  all  of  the  boundary  conditions  are 
implemented  at  the  interface  of  the  guard  cell  and  the  first  cell  in  the  domain.  Therefore,  the 
Dirilecht  boundary  condition  is  implemented  by  setting  the  values  of  the  guard  cells  as  a 
reflection  of  the  internal  cell  values,  6P(nx+l,j)=  5Pg  -  -6P  (nxj). 
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5. 


NUMERICAL  TESTS 


The  solid  combustion  problem  involves  both  pressure  gradient  and  pressure  boundary 
conditions.  For  these  boundary  conditions,  which  may  affect  accuracy  and  convergence,  the 
computer  code  is  tested  by  comparing  the  numerical  results  against  analytical  solutions  for  a  test 
problem. 

Steady  temperature  of  the  electric  heater  (Arpact  [14])  can  be  described  by  the  following 
equations  on  a  rectangular  domain  of  size  2Lx2H  as  shown  in  Figures  3a  and  3b. 


y 

t  u-0 


du/dy=0 


u=0 


(0,0) 

du/dy=0 


Figure  3a  Heat  Transfer  from  a  Square 
Heater 


Figure  3b  Computational  Domain  for 
the  Test  Problem 


The  temperature  field  is  described  by  the  following  equation: 

+  g^u  +  1  =  0.  (5.1) 

dx^  dy^ 

The  boundary  conditions  are 

^  (0,  y)  =  0,  ^  (x,0)=  0,  u(  L,y)  =  0,  u(  x,  H )  =  0.  (5.2a,b,c,d) 

dx  dy 

Equation  (5.1)  is  discretized  on  a  100x100  cell  uniform  grid.  The  resulting  equations  (4.17)  are 
solved  by  Bi-CGSTAB  for  the  temperature  field,  u,  and  are  shown  in  Table  2. 
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Analytical  Numerical  (100x1 00  cells),  tolerance=1 . E-3 


■i 

Case  1 

Case  2 

X 

WM 

u(x,y) 

u(x,y) 

Error 

Error 

%error 

0.1 

0.1 

0.28972 

0.28913 

5.90E-04 

2.30E-03 

0.79 

0.3 

0.1 

0.26958 

0.26852 

1.06E-03 

0.39 

0.27143 

1.85E-03 

0.69 

0.5 

0.1 

0.22754 

0.22598 

1.56E-03 

0.69 

0.22893 

1.39E-03 

0.61 

0.7 

0.1 

0.16012 

0.15796 

2.16E-03 

1.35 

0.16099 

8.70E-04 

0.54 

0.9 

0.1 

0.06227 

0.05935 

2.92E-03 

4.69 

0.06256 

2.90E-04 

0.47 

0.1 

0.3 

0.26957 

0.26852 

1.05E-03 

0.39 

0.27143 

1.86E-03 

0.69 

0.3 

0.3 

0.25118 

0.24974 

1.44E-03 

0.57 

0.25264 

1.46E-03 

0.58 

0.5 

0.3 

0.21264 

0.21080 

1.84E-03 

0.87 

0.21370 

1.06E-03 

0.50 

0.7 

0.3 

0.15035 

0.14806 

2.29E-03 

1.52 

0.15100 

6.50E-04 

0.43 

0.9 

0.3 

0.05888 

0.05601 

2.87E-03 

4.87 

0.05907 

1.90E-04 

0.32 

0.1 

0.5 

0.22753 

0.22598 

1.55E-03 

0.68 

0.22893 

1.40E-03 

0.62 

0.3 

0.5 

0.21263 

0.21080 

1.83E-03 

0.86 

0.21370 

1.07E-03 

0.50 

0.5 

0.5 

0.18116 

0.17909 

2.07E-03 

1.14 

0.18192 

7.60E-04 

0.42 

0.7 

0.5 

0.12949 

0.12719 

2.30E-03 

0.12993 

4.40E-04 

0.34 

0.9 

0.5 

0.05149 

0.04891 

2.58E-03 

0.05165 

1.60E-04 

0.31 

0.1 

0.7 

0.16011 

0.15796 

2.15E-03 

1.34 

0.16099 

8.80E-04 

0.55 

0.3 

0.7 

0.15034 

0.14806 

2.28E-03 

1.52 

0.15100 

6.60E-04 

0.44 

0.5 

0.7 

0.12948 

0.12719 

2.29E-03 

1.77 

0.12993 

4.50E-04 

0.35 

0.7 

0.7 

0.09444 

0.09222 

2.22E-03 

2.35 

0.09468 

2.40E-04 

0.25 

0.9 

0.7 

0.03892 

0.03677 

2.15E-03 

5.52 

0.03897 

5.00E-05 

0.13 

0.1 

0.9 

0.06226 

0.05935 

2.91  E-03 

4.67 

0.06256 

3.00E-04 

0.48 

0.3 

0.9 

0.05885 

0.05601 

2.84E-03 

4.83 

0.05907 

2.20E-04 

0.37 

0.5 

0.9 

0.0515 

0.04891 

2.59E-03 

5.03 

0.05165 

1.50E-04 

0.29 

0.7 

0.9 

0.0389 

0.03677 

2.13E-03 

5.48 

0.03897 

7.00E-05 

0.18 

0.9 

0.9 

0.01754 

0.01608 

1.46E-03 

8.32 

0.01751 

-3.00E-05 

0.17 

Numerical  treatments  for  zero  boundary  conditions: 

Case  1 :  u(nx+1  ,j)  =  -u(nx.j)  Case  2:  u(nx+1  ,j)  =  0. 

u(i,ny+1)=  -u(i,ny)  u(i,  ny+1 )  =  0. 


Table  2  Comparision  of  Bi-CGSTAB  Computations  with  the  Exact 
Analytical  Solution  for  the  Square  Heater 
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Also  shown  in  Table  2  is  the  analytical  solution,  which  is  represented  as  a  series  (n=6)  given  by 

CO 

u(x,y)/L^=0.5*[l-(x/L)^]-22:  (-!)"( cosh  Xn  v  ^  cosOl.  x;>  (5.3) 

(Xn  Lf  cosh  h 

where  ^nL=(2n+l)7i/2,  n=l,2,..6,  and  n  represents  the  number  of  terms  in  the  series.  The 
numerical  solution  was  obtained  for  two  different  boundary  conditions.  Case  1  uses  the 
reflective  boundary  conditions.  In  Case  2  the  guard  cell  values  are  directly  set  to  zero.  Clearly, 
the  numerical  solution  is  generally  within  few  percent  of  the  analytical  solution  with  the 
reflective  boundary  conditions.  Solutions  for  case  2,  however,  have  almost  a  factor  often  less 
error  near  the  domain  comers  than  the  those  for  case  1 .  This  suggests  that  the  exact 
implementation  of  the  boundary  conditions  may  affect  the  accuracy  of  the  solution  near  the 
comers  of  the  domain. 


6.  RESULTS  AND  DISCUSSION 

For  PMMA  combustion,  the  rectangular  domain  shown  in  Figure  1  is  descritized  by 
192x144  finite  volume  cells.  The  cells  are  packed  closely  near  the  leading  edge  region  of  the 
PMMA  plate  and  are  stretched  with  distance  from  the  leading  edge.  The  smallest  cells  are  0.2 
mm  size.  The  gases  are  ignited  for  the  first  2000  time  steps  (time  step=le-05  sec).  At  each  time 
step,  solution  of  equation  (3.20)  and  the  final  solutions  of  equations  (3.1)  to  (3.6)  are  obtained  by 
time  step  splitting. 

The  boundary  conditions  at  the  PMMA  surface  describe  the  no-slip  condition,  the 
interfacial  transport  of  mass,  momentum  and  energy,  and  pyrolysis  kinetics,  all  of  which  are 
given  elsewhere  [1-4].  The  interfacial  temperature,  concentrations  of  the  species,  and  fuel 
ejection  rate  due  to  pyrolysis  are  unknowns  and  are  determined  as  a  part  of  the  solution.  At  each 
time  step,  a  multi-variable  fixed  point  iterative  (MVFP)  method  is  used  to  generate  new  values  of 
the  interfacial  quantities.  The  new  interfacial  values  are  used  as  boundary  conditions  and 
equations  (3.20),  and  (3.1)  to  (3.6)  are  re-solved.  After  convergence  is  obtained  for  the 
interfacial  quantities,  the  computations  are  performed  for  the  next  time  step.  These  iterations  are 
important  since  they  can  significantly  affect  the  local  burning  rates  of  PMMA.  Computations 
were  performed  on  SGI  (Silicon  Graphics)  work  station  equipped  with  300  MHZ,  RIOOOO,  MIPS 
processor. 

Figure  4  shows  the  solutions  of  equation  (3.20)  for  pressure  using  Bi-CGSTAB.  They 
correspond  to  steady  state  (12000  time  steps  of  5e-05  sec).  They  show  a  small  raise  in  5P  near 
the  leading  edge  of  the  PMMA  plate  (x=l  .7  cm).  Combustion  rates  are  highest  near  this  region, 
where  fresh  air  comes  in  contact  with  the  fuel  vapor  generated  at  the  PMMA  surface.  The 
combustion  rates  decrease  sharply  with  the  distance  from  the  leading  edge  region  and  result  in  a 
sharp  decrease  in  6P  as  exhibited  by  Figure  4.  A  much  higher  raise  in  6P  was  found  during  the 
ignition  period  and  it  falls  steadily  with  time. 
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Distance  along  the  Surface,  X,  cm 


Figure  4  Solutions  of  the  Pressure  Equation  using  Bi-CGSTAB  for 
PMMA  Combustion 

Figures  5  and  6  show  the  axial  and  vertical  velocity  components  obtained  from  equation  (3.16). 
Clearly,  the  axial  velocity  decreases  as  the  gases  approach  the  intense  combustion  region  near 
the  leading  edge  of  the  PMMA  plate  and  accelerate  downstream  with  distance  x.  The  axial 
velocity  is  small  near  the  surface  due  to  no-slip,  reaches  peak  value  just  above  the  flame,  and 
approaches  free  stream  value  (84.0  cm/sec)  far  from  the  surface  in  y-direction.  The  vertical 
component  of  the  velocity  shows  a  large  increase  near  the  leading  edge  where  6P  is  the  largest. 
The  vertical  velocity  of  the  fuel  vapors  due  to  PMMA  pyrolysis  is  much  smaller  that  the  peak 
value  shown  in  Figure  6.  Together  these  results  show  that  combustion  near  the  leading  edge 
region  diverts  the  air  flow  upward  away  from  the  PMMA  surface.  This  may  have  important 
implications  on  how  a  suppression  agent  such  as  water  mist  might  distribute  itself  near  this 
region  of  intense  combustion. 
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Figure  5  Axial  Velocity  Contours  for  PMMA  Combustion  Obtained 
by  using  Bi-CGSTAB 
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Figure  6  Vertical  Velocity  Contours  Obtained  by  using  Bi-CGSTAB 

At  the  PMMA  surface  the  local  regression  rate  is  directly  proportional  to  the  local  heat  feed  back 
from  the  flame  to  the  surface.  The  heat  transport  drives  the  PMMA  pyrolysis  that  turns  the  solid 
into  fuel  vapor.  Solutions  of  the  full  Navier-Stokes  equations  are  used  to  evaluate  the  heat  feed 
back.  The  local  regression  rate  along  the  length  of  the  PMMA  plate  is  shown  in  Figure  7.  The 
regression  rate  is  highest  near  the  leading  edge  and  decreases  sharply  within  a  short  distance 
from  the  leading  edge.  The  sharp  decrease  in  the  regression  rate  is  due  to  the  boundary  layer 
structure  of  the  flame.  The  boundary  layer  thickness  increases  sharply  with  distance  from  the 
leading  edge  and  results  in  a  heat  feed  back  profile  that  decreases  with  the  distance. 
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Figure  7  Regression  Rate  Profile  for  the  Boundary  Layer  Combustion  of  PMMA 

Instead  of  using  Bi-CGSTAB,  computations  were  also  performed  using  the  multi-grid 
method  (MGRID,  with  level  4  and  cycle  2)  for  identical  inputs,  which  were  discussed  above. 
Therefore,  all  of  the  computations  are  identical  except  for  the  solution  of  the  pressure  equation 
(3.20).  The  results  are  summarized  in  the  second  and  third  columns  of  Table  3.  The  leading 
edge  of  the  PMMA  plate  is  located  at  1=41  cell  in  the  computational  domain.  Table  3  shows  that 
the  temperature  (K),  pressure  (dynes/cm^),  and  local  burning  rates  (gm/cm^sec)  are  nearly  the 
same  for  both  methods  near  the  leading  edge  (1=42,  j=l)  at  steady  state  (12001  time  steps).  The 
maximum  flame  temperatures,  which  occur  inside  the  domain  (i=64,j=10),  are  also  identical  for 
both  the  methods.  Table  3  also  shows  the  values  of  pressure  at  the  inlet  (i=0),  where  the  pressure 
gradient  is  set  to  zero  (d6P/dx=0),  and  at  the  outlet  (i=193),  where  the  pressure  is  set  to 
atmospheric  (6P=0).  The  computations  with  Bi-CGSTAB  appear  to  satisfy  the  exit  boundary 
condition  better  than  those  with  MGRID  especially  during  the  transient  (6001  time  step). 
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6001th  stq) 


SOOrstep 


10001th  stqj 


12001tiistep 


(fynes/cm 


3.29 


66.1 


-116 


-146 


-0.213 


0.780 


2.03 


3.29 


57.0 


1.25 


-0.021 


0.189 


0.091 


1.02 


3.29 


82.1 


-11.4 


0.704 


0.846 


0.281 


0.588 


Table  3  Comparision  of  Bi-CGSTAB  with  MGRID  Results  for 
PMMA  Combustion 


Computer  time  (clock  time)  profiling  was  also  performed  for  both  the  methods  to 
compare  their  efficiency.  These  are  shown  as  seconds  per  time  step  for  each  stage  of  the 
algorithm.  The  computation  time  for  the  solution  of  the  pressure  equation  (3.20)  is  shown  as 
“Implicit”  in  Table  3.  Clearly,  solving  the  pressure  equation  is  significantly  faster  (factor  of  2.5) 
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with  MORE).  The  number  of  MVFP  iterations  was  found  to  be  roughly  the  same  for  both  the 
methods. 

Table  3  also  shows  the  computations  with  two  different  tolerance  limits  for  the  Bi- 
CGSTAB  iterations.  Increased  tolerance  affected  the  pressure  slightly  but  did  not  significantly 
improve  the  temperature  and  burning  rates.  The  computation  time,  however,  increased 
significantly  with  the  increased  tolerance  as  one  would  expect. 


7.  CONCLUSIONS 

Simulation  of  solid  combustion  problems  involve  an  elliptic  equation  for  pressure,  which 
must  be  solved  multiple  times  at  each  time  step  in  order  to  obtain  local  burning  rate  profiles 
accurately.  A  robust  conjugate  gradient  method  (Bi-CGSTAB)  has  been  implemented  for  the 
solution  of  PMMA  combustion.  The  conjugate  gradient  method  is  shown  to  provide  superior 
performance  in  computation  time  over  a  multi-grid  method  (MGRID)  without  any  loss  in  the 
accuracy  for  the  temperature,  burning  rate,  and  pressure  with  tolerance  limit  set  at  le-03. 
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